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Abstract 

We obtain the first exact solution for the stationary one-dimensional Kardar- 
Parisi-Zhang equation. A formula for the distribution of the height is given in terms 
of a Fredholm determinant, which is valid for any finite time t. The expression is 
explicit and compact enough so that it can be evaluated numerically. Furthermore, 
by extending the same scheme, we find an exact formula for the stationary two-point 
correlation function. 

The dynamics of a growing surface has long been a subject of intense study because of its 
wide appearance in nature and relevance to material science; it is also of great theoretical 
interest as an important example of nonlinear nonequilibrium phenomena [1,2]. For a 
surface growing only through local interaction, a prototypical equation was proposed by 
Kardar, Parisi and Zhang (KPZ) in 1986 [3]. In this Letter, we focus on its one- dimensional 
version, 

dh(xA) A / dh(x,t)\ 2 d 2 h(x,t) , 



dt 2 \ dx J dx 

Here h(x,t) represents the height of the surface at position x G R and time t > 0. The 
last term rj(x, t) is taken to be the Gaussian white noise with covariance {rj(x, t)rj(x', £')) = 
D5(x - x')S(t-t'). 

This KPZ equation has been well known and has already been studied extensively [4-8], 
but interest in this equation has revived [9]. First, the efforts to have a better understand- 
ing of the equation, which are from various perspectives and have lasted for more than two 
decades, finally culminated in the achievement of an exact solution for its height distribu- 
tion [10,11]; second it has turned out that such fine information as the height distribution 
can be measured experimentally. In [12], Takeuchi and Sano employed a turbulent liquid 
crystal as a material and found a clear agreement with theoretical predictions. 
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In the studies of growing surface, the fluctuations of the height are of primary im- 
portance. In many examples one observes the "kinetic roughening." As time goes on the 
surface becomes rougher and rougher, even if it starts from a completely flat substrate. 
This concerns the transient behaviors which are certainly interesting, but even more nat- 
ural and important is the stationary situation. Most systems eventually reach their own 
steady states after relaxation times and our understanding of the nonequilibrium steady 
state (NESS) is in general so limited that concrete results for a simple system should give 
us valuable information. Many papers on the KPZ equation have addressed its stationary 
properties with various methods but all attempts to treat this case exactly have so far 
failed. 

In this article, we provide the first exact results for the KPZ equation in the stationary 
situation. For the KPZ equation, the stationary state is given by the two-sided Brownian 
motion (BM) with respect to x [13]: h(x, 0) = aB_(—x) for x < 0,aB + (x) for x > 0, 
where a = {2v)-^ 2 \D x l 2 and B± (x) are two independent standard BMs, as seen in Fig. 1. 
Since the system is translationally invariant, one can set without loss of generality the 
initial height at the origin to be zero, i.e., h(0,0) = 0. In the following, we will show the 
explicit formulas for the height distribution Pv6b[h(x,t) > s] and the two-point space-time 
correlation functions of the slope (d x h(x, t)d x h(0, 0)) (Fig. 1). The obtained expressions 
are explicit enough so that they can be evaluated numerically. We would say that the main 
achievements in our paper are that we could draw their figures as given in Figs. 2 and 3 
below. 
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Figure 1: The two-sided BM initial condition and the two-point function 



Let us briefly recall the theoretical developments so far. Early works using renormal- 
ization group ideas and other analytic methods have identified the critical exponents. It 
was established that when t is large, the fluctuations of the height grow like 0(t 1//3 ) and 
the nontrivial equal-time correlations in space are observed on the 0(t 2 ^) scale. About a 
decade ago, it was realized that a few growth models admit a mapping to a combinatorial 
problem, which allowed us to study their height distributions and their universal scaling 
limits. In these models some of the space, time, and height are discrete; there is a natural 
ultraviolet cutoff, which makes their analysis somewhat easier. The most surprising at the 
time was that the universal height distribution is the same as that of the largest eigenvalue 
of Gaussian ensembles from random matrix theory [14,15]. 
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But it was only very recently that these developments finally led to the identification 
of the exact height distributions for the KPZ equation itself [10,11]. It contains valuable 
information about finite time behaviors, which has significant relevance to real physical 
systems because a slow approach to the limiting values for the moments has been observed 
in experiments as well. For more recent developments, see, for instance, [16-18]. 

A most important finding from the exact solution is that the KPZ equation represents 
the universal crossover between the linear (A = 0) Edwards- Wilkinson (EW) [19] and the 
nonlinear KPZ universality classes. This is realized, for instance, in the fact that the KPZ 
equation is obtained as a weakly asymmetric limit of the asymmetric simple exclusion 
process (ASEP). In other words, the KPZ equation describes the behaviors of the ASEP in 
a certain time regime [5] . The ASEP is a stochastic system of many particles, each of which 
performs an asymmetric random walk under the constraint of volume exclusion with other 
particles and is one of the most well-studied nonequilibrium discrete models. The KPZ 
equation also captures appropriately low temperature behaviors of some discrete random 
directed polymer models in [20,21], demonstrating the universality of the equation. 

The first exact solution of the KPZ equation was obtained for the narrow wedge initial 
condition, for which the shape of the surface at finite time is a parabola [10, 11]. On the 
other side from the previous studies, it has been well established that the height distribution 
depends strongly on the geometry of the system and its initial conditions [15,22]. In 
particular it is crucial whether the macroscopic shape of the surface is curved or flat. This 
was confirmed in an experiment as well [9]. Hence, it is obviously important to extend the 
exact results to other cases, which is not at all trivial. The analysis in [10,11] was based 
on a previous exact result of the ASEP [23]. The ASEP has been extensively studied, 
but the studies of its fluctuation properties are fairly involved, and there are not so many 
generalizations obtained so far. 

In the meantime the same expression for the narrow wedge initial condition was red- 
erived by using a replica method [20,24]. The replica method is well known to be powerful 
for studying systems with randomness but is also known to possess a serious difficulty re- 
lated to the analytic continuation for the replica number. Interestingly enough, for the KPZ 
equation, though there appear divergent sums during the computations one finally arrives 
at the correct finite expression. Hence, a replica analysis of the KPZ problem should be 
useful for clearer understanding of replica method in general. Moreover, this approach has 
turned out to be more suited for various generalizations. In fact multipoint distributions 
were studied in [25] and, more recently, the flat initial condition was treated in [26]. In the 
present Letter, we study the stationary case by a nontrivial extension of this replica ap- 
proach. By a set of scalings of space, time, and height, x — > a 2 x, t — > 2uaH, h — > -^h, 
we can and will do set v — |, — A — D — 1. 

Now, we start to explain our solutions. The main object in this Letter is the height 
distribution function defined by 

F v± , t (s) = Prob (h(0, t) - t/24 > 7tS ) (2) 

where 7 t = (t/2) 1//3 . Note that this reflects previous studies: The average is given by 
(h(x,t)} = t/24, and the fluctuation scales like 0(t 1/ ' 3 ) from the KPZ scaling. It has not 
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been clear how one can take the average over the BM initial condition directly, or if it is 
possible at all. The novel strategy here in studying the stationary state is that we first 
consider a generalized initial condition and then retrieve the desired information from the 
results for it. The initial condition we consider is 

M*,o)4*-<-;j-^ x i° n ( 3) 

I B + {x) + v + x, x > 0, 

where v± are the strength of the drifts. Once this generalized case is solved, then the 
stationary case can be accessed by taking the v± — > limit. 

Applying the Cole-Hopf transformation, Z(x,t) = e~ h ^ x,t \ the KPZ equation (1) is 
linearized as 

dZ(x,t) ld 2 Z(x,t) , , ry/ . 

This can be interpreted as a problem of a directed polymer in random media; Z is its 
partition function. We consider the iVth replica partition function (Z N (x, t)) and compute 
their generating function G t (s) defined as 

G t (s) = J2 [ -^ W ^(Z N (0,t))e^. (5) 

N=0 

Using the Feynman path integral representation of Z and taking the average with 
respect to the random potential 77, the replica partition function can be written in a form, 
(Z N (x,t)) = (x\e~ HNt \$). Here (x\ represents the state with all N particles being at the 
position x and |<I>), the initial state. For the KPZ equation, the Hamiltonian H N turns out 
to be that of the delta-function Bose gas with attractive interaction [27], 

1 N d 2 1 N 

The eigenvalues and eigenfunctions can be constructed by using the Bethe ansatz [20,24, 
28,29]. For a set of quasimomenta z/s, the eigenvalue E z is given by E z = | J2f=i z j an( i 
the eigenfunction by (x u ■ ■ ■ , x N \^ z ) = C z J2 P( - Sn sgnP Ui< j<k <N { z P(i) ~ z P{k) 
+isgn(xj — Xk)) exp (i Yla=i z P(i) x i^ > where is the set of permutations with iV elements 
and C z is the normalization constant. For the 5-Bose gas with attractive interaction, the 
quasimomenta Zj (1 < j < N) are in general complex numbers and are aligned in the form 
of "strings" [24]. The completeness of these Bethe states was proved very recently [30]. 
We will compute this moment by expanding as (Z N (x,t)) = J2 Z e~ Ezt {x\^ z ) z \§) . 

Because the Brownian motion is a Gaussian process, one can perform the average over 
the initial distribution (3) and the inner product of (\l/ 2 |$) can also be explicitly calculated 
under the condition v± > . It includes the summation over coming from the Bethe 
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eigenstate, which looks highly nontrivial but we have found that there is a combinatorial 
formula which factorizes it [33]. We obtain 



Um=l(- lZ m +V-- \/2){-lZ* m - V + + 1/2) 

Now we can follow the arguments in [31]. We further modify the expression of the replica 
partition function and write it in terms of a determinant by using the Cauchy's determinant 
formula. After some computation, one arrives at an expression for the generating function. 

For the narrow wedge [10, 11] and the half BM initial condition [31], this generating 
function itself is written as a Fredholm determinant, but this time it is not because of the 
novel factor n/li( t '+ + v - ~ 0- The difficulty can be avoided by considering a further 
generalized initial condition for the KPZ equation in which the initial overall height x 
obeys a certain probability distribution. A similar argument was already used for other 
discrete models [22,32,34]. If we denote the height with this initial distribution by h, we 
have h = h — x, where h is the original height for which h(0, 0) = 0. The random variable 
X is taken to be independent of h and hence the A^th moment of e~ h factorizes as (e~ Nh ) = 
(e~ Nh )(e Nx ). Moreover, the probability distribution functions F(s) = Prob[/i(0,t) > 7^] 
and F(s) = Prob[/i(0, t) > -y t s] are related as F(s) = } ljL F(s), where k is the Laplace 

K wt da' 

transform of the pdf of x, and we assume can be Taylor expanded around £ = 0. 

Here it should be noticed that if e x obeys the inverse gamma distribution with parameter 
v + +V-, i.e., if e~ x obeys the gamma distribution with the same parameter, its A^tli moment 
is given by l/n«li( t '+ + v - ~ 0> which exactly compensates for the extra factor above. 
Hence the generating function G t (s) corresponding to h can be written as a Fredholm 
determinant. 

Following the arguments in [20,25,31], one can recover the information of the distribu- 
tion for h from the generating function G t - Combining this and a fact = r(t> + £)/T(t>) 
for the case where e x is the inverse gamma random variable with parameter v, we find that 
the height distribution for the initial condition (3) is given by 



I - / due e7t(s " v v± ,t{ u ) 



(8) 



T(v + + +-it 1 d/ds) 
Here v v ±,t{ u ) ls expressed as a difference of two Fredholm determinants, 

u v± , t (u) = det (1 - P u {B v t - P T ^P U ) - det (l - P u B T t P u ) , (9) 
where P s represents the projection onto (s, 00), 

P£i(6, 6) = Ai£ ^1, ^, v-, v + ^j Ai r r K v+, v^j , (10) 



(6,6) = j dy x _ l e _ lty &il (ti + V, K v-, v + ^j Ai£ + V, ^ v+, v^j , (11) 



5 



and 



A&aA^^l ^il^L, (12) 

1 b 

where T z represents the contour from — oo to oo which, along the way, passes below the 
pole at z = id/b. 

With this expression it is not hard to do analytic continuation for v± to their negative 
regions and then take the v± — >■ limit to study the stationary situation. Our result for 
the height distribution for the stationary KPZ equation is written as 

Prob [h(0, t) - 1/24 > lt s] = / du lt e^ s - u ^- lt(a - u) ^ t {u)- (13) 

r(i + 7 t d/ds) y_oo 

Here vo t t{u) is obtained from v v± .t{u) by taking v± — > limit as in [34]. Details will be 
given elsewhere [33]. By using this formula, one can numerically evaluate the distribution. 
In Fig. 2, we plot this distribution (13) for 7 t = 1. Once we have the formula for finite 
t, it is not difficult to consider the scaling limit. In our setting, it can be studied simply 
by taking t — > oo. One can confirm that it tends to Fq as given in [34]. This function 
already appeared in the studies of lattice growth models and describes the stationary height 
distribution in the scaling limit [34-36]. 




Figure 2: Stationary height distribution for the KPZ equation for ^ t = 1. The dashed 
curve is Fq. The dots (•) indicate the Monte Carlo simulation for WASEP. We set the 
difference between left and right hopping rates 0.15 with their sum normalized to be unity, 
3951 Monte Carlo steps (corresponding to j t = 1), and 10000 samples. 



It is not hard to generalize our analysis to the height near the origin. From the KPZ 
scaling, the nontrivial correlation is expected on 0(t 2//3 ) scale, and we can see that the 
height distribution at x — 2j^X with X finite is written as the same form as in (13) with its 
kernel slightly modified to include X. This can be readily used to obtain the exact formula 
for the stationary two-point correlations. Let us set C(x,t) = ((h(x,t) — (h(x,t))) 2 ). One 
can check that the second derivative of this function with respect to x is, in fact, twice the 
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two-point correlation: d 2 C{x,t) = 2(d x h(x, t)d x h(0, 0)). From the KPZ scaling, we also 
introduce the scaled version, g t ( y ) = (2t)- 2 / 3 C ((2t) 2 ^y,t). Plotted in Fig. 3 is the second 
derivative of this function gt(y)- In the scaling limit, it tends to the function g(s), which 
was studied before in [34-37]. More discussions will be given in [33]. 

Here we remark that in statistical mechanics, computing the dynamical correlation 
functions is one of the most important objectives of the theory. For instance, in the 
Kubo formula, one needs the dynamical current-current correlation function in equilibrium. 
This is already known to be very difficult in general and it is even more interesting and 
challenging to compute the dynamical correlation for nonequilibrium systems, in particular 
in its stationary state. In addition, the Fourier transform of the two-point function is known 
as the structure function, which can be measured directly using (e.g., neutron) scattering 
experiments. 




Figure 3: Stationary two-point correlation function g' t '(y) for the KPZ equation for j t = 1. 
The dashed curve represents the one in the scaling limit g"(y). The dots indicates the 
Monte Carlo result of the WASEP, where parameters are the same as those in Fig. 2. 

Our exact results represent the universal crossover for the KPZ equation with t being 
a parameter. As explained above, in a certain time regime, the fluctuation properties of 
the weakly ASEP (WASEP) are described by the KPZ equation. We confirmed this by 
comparing our theoretical predictions with Monte Carlo simulations of the WASEP (Figs. 
2 and 3). The dots in Figs. 2 and 3 represent the Monte Carlo results for the current 
distribution and density-density correlation in the WASEP respectively. 

To summarize, we have found the first explicit formulas for the height distribution and 
the two-point correlation functions for the KPZ equation in its stationary regime. They are 
explicit enough to allow numerical evaluations as depicted in Figs. 2 and 3. One observes 
nontrivial deviations from the scaling limit. They contain valuable information about the 
finite t corrections which should be observed in real experiments. One possibility to realize 
the stationary situation in an experiment is to use again the liquid crystal turbulence 
because the system is fairy flexible for setting a desired initial condition by adjusting a 
shape of the laser pulse. It might also be possible to observe our prediction by performing 
a careful experiment of paper combustion for which the stationary situation was already 
studied in a previous publication [38]. 
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